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Abstract 

We discuss the hard-hexagon and hard-square problems, as well as the correspond- 
ing problem on the honeycomb lattice. The case when the activity is unity is of interest 
to combinatorialists, being the problem of counting binary matrices with no two adja- 
cent l's. For this case we use the powerful corner transfer matrix method to numerically 
evaluate the partition function per site, density and some near-neighbour correlations 
to high accuracy. In particular for the square lattice we obtain the partition function 
per site to 43 decimal places. 

1 Introduction 

Partition function of hard squares 

In recent years there has been interest amongst combinatorialists in the problem of 
counting the number of "legal matrices" [||, 10, 11, |l2| . These are matrices whose 
elements are either or 1. Two elements are said to be 'adjacent' if they lie in positions 
and (i + 1, j), or if they lie in positions and + 1). No two l's are allowed 
to be adjacent. 

In statistical mechanics this is known as the 'hard squares 'problem: how many ways 
are there of putting particles on a square lattice of N sites so that no two share the 
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same site, or are adjacent. (An economical way of formulating this rule is to say that 
there is at most one particle on every edge.) 

In statistical mechanics we usually generalize this problem by specifying the number 
of particles. Thus we might define g(m, N) to be the number of ways of putting m 
particles on a given lattice of N sites, subject to the above exclusion rule. This is the 
'canonical' partition function. The 'grand canonical partition function' is 

E N (z) = 5>(m,A0z m , (1) 

m 

z being an arbitrary real (or complex) variable. For z real and positive we expect on 
physical grounds that the limit 

k(z) = lim [E N (z)) l / N (2) 

will exist and be independent of the shape of the lattice, so long as the shape is not 
pathological - it must certainly become infinite in all directions. 

It follows that what combinatorialists are interested in is H_/v(l) and /-t(l). While the 
function k{z) has been extensively investigated, by series expansions and numerically 
J5], 15, 25, |26|, less attention has been paid by statistical mechanics to its value at 



z 



1. Estimates for k(1) of 1.5030 and 1.503048082 are given in |2§ and |2j], and 



lower and upper bounds of 1.503263 and 1.517513 in [21]. 

Amongst combinatorialists, Engel |l^] obtained the estimate 1.50304808, and Calkin 
and Wilf || have recently obtained 1.5030480824753323, with rigorous lower and up- 
per bounds 1.503047782 and 1.5035148. This estimate has been refined by McKay to 



1.50304808247533226 |2] 



Here we use the corner transfer matrix method to calculate k(1) and various related 
quantities and obtain the result 1.5030480824753322643220663294755536893857810. 
Based on the observed convergence of the sequence of successive approximations (and 
on the expected error of A.%), we believe this result to be correct to the 43 decimal 
places given. 
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Related models and quantities 

For the corresponding problem on the triangular lattice, known as 'hard hexagons', 
Metcalf and Yang in 1978 (2^] conjectured that log/-c(l) = 1/3. This intriguing conjec- 
ture prompted Baxter and Tsang ||] to obtain more accurate numerical values, using 
the corner transfer matrix method ||. They disproved Metcalf and Yang's conjecture, 
but as a result the author stumbled on the fact that the problem could be solved exactly 
by the Yang-Baxter' method (C.N. Yang, rather than his brother CP.). |2|, || 

For the square and honeycomb lattices, the problem has not been solved exactly 
and the suggestive clues of hard hexagons are missing. Even so, the method provides 
a powerful way of calculating k(1) to considerable numerical precision. We do this 
here, obtaining it to 42 and 39 decimal places for the square and honeycomb lattices, 
respectively. For completeness, we also evaluate k(1) for the triangular lattice, using 
the known exact result, and give it here to 55 decimal places. 

A related quantity is the density, or probability that a given site is occupied: 

p(z) = E N (zr 1 J2HN)g(m,N)z m = z^ln K (z) , (3) 

m aZ 

also expected to tend to a limit for a large lattice. This is a 'single-site' property and is 
also readily calculated by the corner transfer matrix method, so we present numerical 
values for p(l) here, along with some next-nearest neighbour correlation results. 

It should be noted that the models are expected to undergo a phase transition from a 
disordered fluid phase to an ordered solid crystalline phase at z = 11.09017..., 3.7962..., 
7.92..., p = 0.27639..., 0.368..., 0.422..., for the triangular, square and honeycomb lat- 
tices, respectively, p], 0, |2?j]. For each lattice, these critical values of z are considerably 
greater than z = 1, so the gas is well and truly in the fluid phase and we expect suc- 
cessive finite truncations of the corner transfer matrix equations (with no sub-lattice 
symmetry breaking) to converge rapidly, and indeed we observe this. The method 
converges more slowly as z approaches the critical value. 
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2 Low density series expansions 

The series expansions were obtained for these models in the sixties |l7], [l8|, |28| to orders 
z 8 , z 13 , z 13 , respectively. They have since been extended, e.g. in j|, p. 408], @. Here 
we give each to order z 13 : 

triangular : k = 1 + z - 3 z 2 + 16 z 3 - 106 z 4 + 789 z 5 - 6318 z 6 + 

53198 z 7 - 464673 z 8 + 4174088 z 9 - 38332500 z 10 + 358371345 z 11 - 
3400238419 z 12 + 32664062740 z 13 + 0{z u ) 

square : « = 1 + z - 2 z 2 + 8 z 3 - 40 z 4 + 225 z 5 - 1362 z 6 + 

8670 z 7 - 57253 z 8 + 388802 z 9 - 2699202 z 10 + 19076006 z 11 - (4) 
136815282 z 12 + 993465248 z 13 + 0(z 14 ) 

honeycomb: k 2 = 1 + 2 z - 2 z 2 + 6 z 3 - 23 z 4 + 100 z 5 - 469 z 6 + 

2314 z 7 - 11841 z 8 + 62286 z 9 - 334804 z 10 + 1831358 z 11 - 
10162679 z 12 + 57080840 z 13 + 0(z 14 ) 

3 Corner transfer matrix equations 

Label the sites of the lattice i = 1, . . . , N. With site i associate an occupation number 
(or spin) <jj, with value if the site is empty, 1 if it is full. Then for all three lattices 
the partition function is 

^n(z) = En^ 1 n <wi) > ( 5 ) 

i <ij> 
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where the sum is over all values or 1 of all the occupation numbers (without restric- 
tion), the second product is over all edges < i,j > of the lattice, and 

e(0,0) = e(0,l) = e(l,0) = 1 , e(l,l)=0 . (6) 
For the triangular lattice the corner transfer matrix equations are given in || . They 

are 

J2F(a,b)A 2 (b)F(b,a) = ^A\a) , 

b 

J2w{a,b,c)F(a,c)A(c)F(c,b) = rjA(a)F(a,b)A(b) , (7) 

c 

k = r, 2 /i . 

Here a,b,c are site occupation numbers, taking the values or 1, and w(a,b,c) is the 
weight function of a triangular face of the lattice, which for our hard molecule model 
is 

w(a,b,c) = z^ a+b+c ^ 6 e(a,b)e(b,c)e(c,a) . 

For the square lattice the corner transfer matrix equations are given in jj| and ||]. 
Specializing to the translation-invariant case (with no symmetry breaking), they are 

J2F(a,b)A 2 (b)F(b,a) = U\a) , 
b 

J2w(a,c,d,b)F(a,c)A(c)F(c,d)A(d)F(d,b) = r 1 A(a)F(a,b)A(b) , (8) 

c,d 

K = r]/i , 

the face weight function being 

w(a,b,c,d) = z (a+b+c+d)/4 e{a,b)e(b,c)e(c,d)e(d,a) . 

For the honeycomb lattice, we found it convenient to transform the model to one 
on the triangular lattice by dividing it into two sub-lattices and summing over all 
the spins on one sub-lattice. (This is a decimation transformation.) The resulting 
triangular lattice model has different face weights w\ , Wi for up-pointing and down- 
pointing triangles, respectively: 

iDi(a,6,c) = z( a+6+c )/ 3 + zS a0 5 b o5 c0 , w 2 (a,b,c) = 1 , 
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and the equations (|7|) generalize to 

^Fi(a,b)Mb)M b ) F i( b > a ) = ^A j (a)A i (a)A j (a)A i {0') , 

b 

Y J Wi{a,b,c)F i (a,c)A i {c)F i (c,b) = ^^(0)^(0,6)^(6) , (9) 

c 

« = (W0 1/2 > 

for (i,j) = (1,2) and (2,1). Here k is the partition function per site of the original 
honeycomb lattice. 

In each case A(a) is the corner transfer matrix centred on a site of occupation 
number a. Define the matrix 

M(a)=A 6 (a), A 4 (a) , [A 1 (a)A 2 (a)} 3 , (10) 

for the three lattices, respectively. Then the density is 

p = Trace M(l)/[Trace M(0) + Trace M(l)] . (11) 

Similarly, F(a, b) is the 'half-row' transfer matrix associated with an edge with occu- 
pation numbers a, b. For this reason F(l, 1) is the zero matrix for the first two lattices. 
(For the honeycomb lattice the decimation transformation means that Fi(a,b) is as- 
sociated with next-nearest neghbour sites a, 6, so i^(l, 1) is not zero.) Other local 
occupation probabilities can be obtained by building up a face or faces of the lattice. 
For instance, for the square lattice the probability that two next-nearest-neighbour (i.e. 
diagonally adjacent) sites are simultaneously occupied is 

w(l, 0, 1, 0)Trace [A(1)F(1, 0)A(0)F(0, 0)^(0)F(0, 1)] 
Pnn ~ Eabcd w{a, b, c, (i)Trace [A(a)F(a, b)A(b)F(b, c)A(c)F(c, d)A(d)F(d, a)} 

(12) 

It can be convenient to define the larger matrices 
/ 



A 




F I FM 1 (13) 

F(1,0) F(l,l) 



M(0) 
M(l) 
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(Generalizing as needed for the honeycomb lattice.) 

These extended matrices A, A4 are real and symmetric. The corner transfer 
matrix equations above define these matrices to within irrelevant similarity transfor- 
mations: it is convenient to remove this arbitrariness by choosing A and A4 to be 
diagonal. 

The equations are exact in the limit when the matrices are infinite-dimensional. 
Finite-dimensional truncations provide a sequence of rapidly converging approxima- 
tions and the techniques needed to handle these are explained in [jl], ^, g|. For each 
lattice a column k, b of the second equation plays the role of an eigenvalue equation, 
the eigenvalue being the diagonal element k, b of A, and the eigenvector being column 
k, b of T ' . The first equation fixes the normalization of the eigenvector. In any finite 
truncation there are more eigenvalues than needed: one selects those giving the larger 
eigenvalues of A4. One can normalise the diagonal matrix hA to have maximum el- 
ement unity, the diagonal elements being arranged in numerically decreasing order. 
These elements do not change dramatically as one increases dimensionality, so one can 
regard the truncation where A4 has dimensionality r as providing an approximation 
for the largest r elements of the true infinite-dimensional matrix A4. The magnitude 
of the largest element omitted provides a measure (in practice it seems a somewhat 
over-confident measure) of the accuracy of any given truncation. 

For each lattice, and for any finite truncation, the CTM equations can be obtained 
from a variational principle for k. 

4 Simple approximations 
Bethe approximation 

A simple approximation for lattice models is to replace the lattice by a Bethe lattice 
(i.e. a Cayley tree with surface effects excluded) of the same coordination number q. 
Specializing the formula given in Chapter 4 of we get 

k = (l + m)/(l -m 2 ) q/2 , p = m/(l+m) , (14) 
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where m is related to z by 

z = m/{l-m) q . (15) 

This approximation is self-consistent in that the relation (0) is satisfied. When z is 
small, so is m and we readily obtain the series expansions for the three lattices (with 
q = 6,4,3, respectively): 

k = 1 + z - 3z 2 + 18 z 3 + 0(z 4 ) , 

k = l + z-2z 2 + 8z 3 -41z 4 + 0(z 5 ) , (16) 
k 2 = l + 2z-2z 2 + 6z 3 -23z 4 + 100z 5 - 470z 6 + O(z 7 ) . 

Let q' be the number of edges round a face of the lattice (this is the coordination 
number of the dual lattice), so q' = 3, 4, 6 for the three lattices, respectively. Comparing 
( |l~6| ) with (|j), we see that for each lattice the Bethe approximation fails at order z q ' , 
which is the order at which a circuit first makes a contribution to the partition function. 



Kramers- Wannier approximation 

A more accurate approximation, which at least in some circumstances is related to 



that of Kramers and Wannier [20], is to take the corner transfer matrices A(a),F(a, b) 
in ([?]) - (H) to be one-by-one matrices. Then A, T are two-by-two matrices. For the 
triangular lattice this gives 

K = (1 - m) 3 /(l - 2m) 2 , p = m/(l+m) , (17) 

where Ai has the diagonal entries 1, m and 

z = m{l - mf/(l - 2mf . 

Expanding, this gives 

K = 1 + z - 3 z 2 + 16 z 3 - 106 z 4 + 789 z 5 - 6319 z 6 + 0{z 7 ) . (18) 
For the square lattice, 

K = (l + t) 2 /{l + t-t 2 ) , p = m/(l + m) , (19) 



where now t, m, z are related by 

z = t(l + t) A /(l+t-t 2 ) , m = t/(l+t-t 2 ) , 
and t is small when z is. This yields the expansion 

k = 1 + z - 2 z 2 + 8 z 3 - 40 z 4 + 225 z 5 - 1362 z 6 + 8670 z 7 - (20) 
57254 z 8 + 388830 z 9 + 0(z w ) . 

For the honeycomb lattice, truncating Ai{a) and Fi(a,b) in (|9|) to one- by-one, we 
obtain eight distinct equations for £,771,772 an d the elements of the matrices. We have 
not succeeded in significantly simplifying the solution of these eqns, but we do note 
that to leading order in z we can naturally choose, 

At = 





and 

( 1 zV3 

Ti = T 2 = \ 

\ z 1 ' 3 z 2 / 3 

Without loss of generality we can fix Ai to be the unit matrix, and normalise A2, T\ , J~2 
to have top-left element equal to unity. We can then iteratively solve the eight equations 
for the eight remaining unknowns to successively higher orders, obtaining: 

k 2 = 1 + 2 z - 2 z 2 + 6 z 3 - 23 z A + 100 z 5 - 469 z e + 

2314 z 7 - 11841 z 8 + 62286 z 9 - 334804 z 10 + 1831358 z 11 - (21) 
10162680 z 12 + 57080872 z 13 + 0{z 1A ) 



We note that for each lattice this 'Kramers-Wannier' approximation is more accu- 
rate than the Bethe approximation, being correct for orders up to (but not including) 

7 2q' 
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One obtains successively more accurate approximations by taking the matrices A, F 
to be larger. This fact is used in |J and W to obtain series expansions. Here we use 
it to obtain accurate numerical values of k and p for z = 1. If ^4(0) is of size no by no 
and ^4(1) is n\ by m, then F(0, 1) is no by rt\ and A, T are n by n, where 

n = no + n\ . (22) 

From now on we re-arrange the combined matrices A, J 7 , S, Ai so that the entries of 
the diagonal matrix A (A1A2 for the honeycomb lattice) are in numerically decreasing 
order. Then S is a diagonal matrix with entries or 1, corresponding to whether the 



centre site is empty or full for that state. The formula (11) can still be written 



Trace S M /Trace M . (23) 



5 Values of k and p for z = 1 

Throughout this section we take z to be 1. For all three lattices the 'Kramers- Wannier' 
approximation has no,ni,n = 1,1,2. The next approximation has no,ni,n = 2,1,3. 
We present the numerical results obtained for the Bethe approximation and these two 
simple corner matrix approximations, together with those for no,ni,n = 3,2,5. It is 
apparent that the results are converging rapidly towards a limit. For the square lattice 
we have continued the sequence of approximations to n = 48, i.e. we have kept the 
largest 48 diagonal elements of A. Then no,ni = 29,19. (In fact, since writing this 
paper we have extended the square lattice calculation to n = 60, working to 55-digit 
accuracy, which has increased our confidence in the numerical results for this case.) 
Similarly for the honeycomb lattice we have continued to n = 20, keeping the largest 
20 elements of A1A2, with no,m = 11,9. 

For the benefit of anyone who wants to repeat these calculations, in the Appendix 
we have given the diagonal elements of A and A1A2 for the square and honeycomb 
lattices with n = 48 and n = 20, respectively. In each case we have included the 
next-largest eigenvalue, which is obtained during the iterative solution procedure, but 
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not actually used in the equations (@) , @ . This provides a measure of the accuracy of 
the approximation. The values of these elements change as one changes n, but not by 
more than one per cent. For instance, the fifth eigenvalue of A for the square lattice 
occurs first when n = 5. For this and the n = 6, 7, 8, 48 truncations it is 0.0077696234, 
0.0077684372, 0.0077696114, 0.0077696231, 0.0077696235. Any given element of A or 
J- tends quite quickly to its limiting n = oo value as n increases, once n is big enough 
for the element to occur at all in the truncation. 

For the triangular lattice the results for n = 2,3,5 are given in Q. We did not 
pursue the corner transfer matrix numerical calculations any further for this lattice 
as (a) the exact result is given in 0], in chapter 14 of 0] and (in algebraic form) 
in |19[ |, (b) the fact that it is exactly solvable is connected with the fact that A has 
degenerate eigenvalues, at any rate in the infinite-matrix-size limit. These degeneracies 
actually complicate the numerical calculation, producing extra degrees of freedom that 
need to be fixed. Unfortunately we observed no such degeneracies for the square and 
honeycomb lattices, so we have no reason to suppose that they are solvable by the 
means used for the triangular case. 

The 'n = oo' results we quote for the triangular lattice are obtained from the 
formulae (14.1.20) and (14.5.14) of @, using the solution of (14.1.18) to 55 significant 
digits: 

x = -0.2549635631051309947933138248965459184034888000019327401 . (24) 

The solutions of ( |i~4|) for the triangular, square and honeycomb lattices (with q = 
6,4,3, respectively) are m = 0.2219104,0.2755080,0.3176722. 

The solution of is m = 0.1932944673, that of (|l|) is t = 0.3598273461, while 



for the honeycomb lattice the 'Kramers-Wannier' two- by- two solution of (|9|), (13) is 
f = 1.2736933548, 771 = 2.4537675065, rj 2 = 1.2413585145, 
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1.395217 


1.502928407 
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1.3954858818 


1.5030479990 


1.546440707581 
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1.3954859724543 


1.50304808247401 


1.54644070878756097 




oo 


1.3954859724791398 


1.5030480824753323 


1.546440708787561419 



Table 1: Values of k for the three planar lattices in the Bethe and Kramers- Wannier (KW) approximations, 
and using higher n by n corner transfer matrix truncations. The extrapolated values (exact to the accuracy 
given) are in the last row. 



1 0.6736226737 
0.6736226737 0.4800087072 



1 0.5940391501 
0.5940391501 0.5960317063 



The values of k and p obtained are given in Tables 1 and 2, respectively. They are 
given to sufficient accuracy to see the convergence of the successive approximations. The 
deviation of an approximation from the extrapolated limit of the sequence (or what is 
almost the same: the deviation from the next approximation) is indeed found to be of 
the order of (or at most two or three orders greater than) the magnitude of the largest 
diagonal element of M. omitted. This gives us some confidence that the sequence is 
indeed converging to the correct value, and the accuracy of the best approximation. 

In fact we have obtained k and p to considerably more accuracy than we are able 
to display in Tables 1 and 2. For the triangular, square and honeycomb lattices they 
are: 



k = 1.3954859724793027352295006635668880689541037281446611908 (55) 
= 1.5030480824753322643220663294755536893857810 (43) 
= 1.54644070878756141848902270530472278026 (38) 
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0.22657081546271469 
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Table 2: The values of p for the three lattices in the various approximations. 



p = 0.1624329213974881529255929066818976201010622684352448332 (55) 
= 0.22657081546271468894199226347129902640080 (41) 
= 0.2424079763616482188205896378263422541 (37) 

The bracketted numbers are the number of decimal places given: for the triangular 
and square lattices we believe the results are accurate to this number of digits. For the 
honeycomb lattice the last two or three digits should be treated with caution. 

Next-nearest neighbour correlation 

For all these models there is if course zero probability that two adjacent sites are 
simultaneously occupied, i.e (<Ji<Jj) = if sites i,j are adjacent. For the square and 
honeycomb lattices, our corner transfer matrix approximations also immediately give 
the probability that two next-nearest-neighbour sites are simultaneously occupied. Yet 
higher correlations can be obtained by constructing the partition function of a small 
lattice, using the A and F matrices appropriately to include the contributions of the 
boundary sites and edges. For the square lattice, with k,m being the sites shown 
in Fig.l, we find (to 38 decimal places) 
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Figure 1: Four sites k, m of the square lattice. 



p 2 = {(Ticrj) = 0.08198200258221599221104301382489246711 

P2 = ((JiCTk) = 0.06967964946634484029326972372242595563 

p 3 = (aicrjcrk) = 0.03024840234446871928473381972162222576 

p 4 = (piOjG k G m ) = 0.01313119289260936136017735696986128175 . 

This information is sufficient to calculate all correlations between the four sites i, j, k, m. 
Since the centre site p is connected to the rest of the lattice only via k,m, we can 
use a "star-square" relation to calculate (a p ) from these correlations [14, eqn. 80], 
[13, |l^, |i~q| . The result is the relation 

6p - 4p 2 - 2 P2 ' + 4p 3 - Pi = 1 • (25) 

This is a useful check on our numerics: it is indeed satisfied to the available accuracy. 

For the honeycomb lattice, with i,j,k,m being the sites shown in Fig. 2, to 33 
decimal places we obtain 

p 2 = (a-io-j) = 0.079618322060174417313883811021777 
p 3 = (oia-jcrk) = 0.026815084372282157838703243933620 
p 3 ' = (aiajam) = 0.032467404849412503526876270338554 . (26) 
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Figure 2: Four sites i, j, k, m of the honeycomb lattice. 

The site p is connected to the rest of the lattice only via i, j, k, so this time we can use 
the "star-triangle" relation to obtain 

5p-3p 2 + P3 = 1 , (27) 
which is indeed satisfied to the available accuracy. 
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Table 3: Diagonal elements of A for the square lattice for the n — 48 approximation 
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Table 4: Diagonal elements of A1A2 for the honeycomb lattice for the n — 20 approximation 
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